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We study the quantum melting of the two-dimensional Wigner crystal using a fixed node quantum 
Monte-Carlo approach. In addition to the two already known phases (Fermi liquid at large density 
and Wigner crystal at low density) , we find a third stable phase at intermediate values of the density. 
The third phase has hybrid behaviors in between a liquid and a solid. This hybrid phase has the 
nodal structure of a Slater determinant constructed out of the bands of a triangular lattice. 



The physics of a system of N electrons confined on 
a two dimensional surface S' is a textbook problem at 
the root of a very large body of literature. Two com- 
peting energies, electrostatic and kinetic, give rise to a 
rich phase diagram. The physics is controlled by the di- 



mensionless parameter r^ 



* e^ I {hr i.^/wn) which is the 



ratio of the average distance between electrons over the 
effective Bohr radius (e is the electronic charge, e the 
dielectric constant, m* the effective mass and n = N jS 
the electronic density). At large density (low r^), the 
kinetic energy dominates and the system is in a Fermi 
liquid phase [jj]. Since the work of Wigner j^] in 1934, it 
is also known that at low density (large r^), the Coulomb 
repulsion dominates and the electrons crystallize onto 
a (Wigner) triangular crystal P, |3]- In their pioneer- 
ing work in 1989, Tanatar and Ceperley ^ were able 
to locate that the quantum melting of the crystal oc- 
curs for a critical value of rs k, 37 ± 5. There work, 
which used a Fixed Node Quantum Monte-Carlo |a| (FN- 
QMC) technique, was followed by more precise numer- 
ics [g and a better description of the liquid phase [3, Q 
that included backflow corrections. 

This simple picture of a, presumably first order, direct 
transition between the solid and the liquid phase is to 
be contrasted with other aspects of the physics of the 
Wigner crystal which show more complex behaviors. For 
instance, its magnetism is believed to include a spin liq- 
uid phase in addition to the ferromagnetic phase found at 
very large r^ 9] . The fcrmionic statistics of the electrons 
is also known to play a crucial role for Ts < 60 where 
the melting of the bosonic Wi gner crystal occurs [ill . 
Also, the classical melting [ll|, |l2, llJ, |lj| (as a func- 
tion of temperature) occurs in two steps. The system 
first looses its translational order but retains some orien- 
tational order (hexatic phase [lj|) while at higher tem- 
perature, all order disappears. The possibility that the 
melting of a quantum crystal would also take place in two 
steps, leading to a highly correlated intermediate phase 
has been discussed as early as 1969 by Andreev and Lif- 
shitz [13, who proposed that a hquid of defects would 
exist together with the crystal state. This proposal has 
been revisited recently in small systems using exact di- 
agonalization techniques [la, [l3| as well as in screened 
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FIG. 1: Energy difference iSuquid — iScrystai (circles) and 
-Ehybrid — -Ecrystai (triaugles) as a function of ra for 72 elec- 
trons in a 48 X 84 grid. Inset: Energies of the three phases at 
v — 1/56 and U — 20 (rs ~ 42.2) as a function of the number 
of particles up to 200 electrons in a 80 x 140 grid. 



systems [l8J . 

In this letter, we study a new phase which is a hybrid 
of a liquid and a solid using a FN-QMC technique sim- 
ilar to the one used in [J|. The FN-QMC approach is 
a very powerful tool to tackle this problem, but it is of 
primary importance to understand the nature of the ap- 
proximations which it involves. The method lies halfway 
between a black box and a variational approach. Techni- 
cally, the FN-QMC algorithm is fed with a wave-function 
called the guiding wave function (GWF) that has to be 
given explicitly, and that should be close to the ground 
state of the system. The FN-QMC algorithm modifies 
the GWF to become as close as possible to the ground 
state of the system, given the constraint that the sign of 
the wave function remains unchanged at every point of 
the Hilbert space. The method gives the best wave func- 
tion for a given structure of the nodes of the GWF and is 
in this sense variational [l9j . Our main result is summa- 
rized in the stability diagram Fig. ^ where the energies 
of the different phases (i.e. associated with the different 
GWFs) are plotted as a function of r^ . The hybrid phase 
is found to be stable in the (critical) region r* < rs < r** 



with r* fs 30 and r** « 80. 

Model. We consider a system of N spinless electrons on 
a square L^ x Ly grid with periodic boundary conditions 
whose Hamiltonian is given by, 



{f,f') 



-E 

2 f-: 



V{f-r')npnp+X. (1) 



The operator cl {cp) creates (destroys) an electron on 
point r with the standard anticommutation relation 
rules. The sum X](r rO ^^ done on the nearest neighbor 
points on the grid and t is the corresponding hopping 
amplitude. The density operator reads np — clcp. U is 
the effective strength of the interaction. The two body 
interaction V{r) is obtained from the bare Coulomb in- 
teraction using the Ewald summation techniques to avoid 
finite size effects, and reads 



y(f) = ^. 



-Erfc(fcc 



L\ 



(2) 



+ 



^ ^Erfc(|i?|/(2/cc)) cos(i? • r). 



2-K 






In the previous equation, kc is a (irrelevant) cut off. The 
vector L takes discrete values L = {nxLx,nyLy) with Ux 
and Uy integer numbers. The vector K also takes discrete 
values, K = {^nx.^riy) and {nx,ny) ^ (0,0). The 

complementary error function is Erfc(r) = -7= J^ e~* dt. 
In order to assure electrostatic neutrality we add a pos- 
itive continuous background (the positive background 
charges are not put on the grid but lie in the continuum). 
The constant term A hence reads, 

X/N ^4t + UV{6) - Ui^^ - C/^ (3) 

kc VTT 



where u 



N 



is the average electronic density. All the 



energies in the problem are measured in unit of N2'Kvt. 
The Tg parameter for this model reads, r^ = U l(2t^f¥u). 
When V ^ \ the role of the grid becomes irrelevant 
and Eqjpi tends toward the continuous model studied 
in Ref fj| provided our energies are multiplied by 2/r^. 
As we shall see however, the presence of the grid gives the 
possibility of constructing new types of GWF. In our nu- 
merics we have used v = 1/56 and v = 1/780. Standard 
two dimensional gas in GaAs heterostructures where the 
underlying grid is given by the Ga and As atoms cor- 
respond to i^ « 1/1000 or 1/ w 1/10000 for the most 
diluted ones. In order for the Wigner crystal to fit into 
the system without distortion, we chose Ly « voLj, and 
N — 2P^ with P integer. The electrons in our study are 
fully spin polarized which corresponds to a system with a 
strong in plane magnetic field. However our results also 



extend to zero field systems since at r^ > 20 the polar- 
ized fluid is more stable than the non polarized one |6|,l31. 
Last, we have added a very small (irrelevant) disorder to 
the system in order to lift the degeneracies of the non 
interacting problem. 

FN-QMC Method. The operator e"^* is apphed 
stochastically to an initial GWF in order to project it 
to the exact ground state. Our implementation is based 
on the Green Function Monte Carlo for lattice Hamil- 
tonians introduced in |20|. Important sampling J21| and 
Fixed Node are implemented as in [ijj by replacing H by 
an effective Hamiltonian i/pN that depends on the GWF. 
if FN forbids the sign of the wave- function to change. The 
energies calculated with iipN are larger than the one of 
the true ground state but smaller than the variational en- 
ergy associated with the guiding wave- function [l9(. At 
V <^1 the technique is equivalent to the continuous fixed 
node diffusive Monte-Carlo used in |j| . The algorithm to 
update the Slater determinants can be found in |23| . By 
sampling directly the time spent by the walkers at one 
point of the Hilbert space using the algorithm described 
in [2OJ we can use arbitrary small time steps and effec- 
tively work in continuous (imaginary) time. Instead of 
using branching, the control of the walkers population is 
done using a fixed number of walkers and the reconfig- 
uration algorithm introduced by Sorella |23|. This algo- 
rithm allows to avoid the bias introduced in the branch- 
ing technique by artificially controlling the walker pop- 
ulation. Quantum averages of physical quantities (...) 
are calculated using the forward walking technique |23l |. 
and hence do not suffer from the bias of mixed estimates. 
A typical point for 72 particles involves 20 independent 
Monte-Carlo runs with 5000 walkers each. 

Guiding wave functions. The GWFs used in our cal- 
culations are Slater determinants multiplied by Jastrow 
functions, 

*(ri,r2...rW) = Det [0,(fj)] x [] J(|7-, - r^l). (4) 

The Jastrow part takes Coulomb interaction into ac- 
count by introducing correlations between electrons. It 
has no nodes, and thus is irrelevant in the FN-QMC re- 
sults. We use modified Yukawa functions J2J|, J(r) = 
—^^(1 — e^'^'^''"'^''), where the distances are measured in 
unit of the average distance between nearest particles and 
A{rg) and B{rg) are (optimized) variational parameters. 
We checked that the FN-QMC results are not sensitive 
to the choice of the Jastrow function. The Slater deter- 
minant of one-body wave functions, Det [(j)i{rj)\ enforces 
the antisymmetric nature of the fermionic wave function 
and is responsible for the nodal structure of the GWF. 
The GWF used in the hterature are constructed out of 
plane waves 4>i{rj) oc e'^''^'^^ for the liquid GWF ^uq and 
localized orbitals 4>i{rj) oc e"''"^""*-' /'^o for the crystal 
GWF *ciy Here the Ui with i ^ {I. ..N} stand for the 



positions of the electrons in the classical crystal and do is 
a variational parameter, ^nq (^cry) provides the exact 
ground state of H at very large (low) density. 

Hybrid GWF. Below we give the detailed construction 
of a new GWF, \l'hyb, such that the (pii^j) are the Bloch 
states of a triangular crystal. First, an effective one-body 
Hamiltonian H^s is constructed for an effective hole in a 
periodic potential given by a classical Wigner crystal. 



-f^cff ^~tJ2 4c?' -U*J2 Wir)np 



(5) 



{r^r') 



where the one-body potential is W{'r) — J2i=i V{r — Ui). 
The singularity of W{r) at r ^ Ui has been removed by 
setting W{ui) = W{ui + (1, 0)) and we checked that our 
results are unaffected by this choice. In a second step, we 
take advantage of the presence of the underlying grid and 
HeS is numerically diagonalized using Lanczos algorithm. 
The A'' orbitals of lowest energy 4>i{r) {1 < i < N) are 
then used to construct the Slater determinant. U* is a 
variational parameter. 

The underlying idea behind the construction of ^hyb is 
to put on the same level the melting of the Wigner crystal 
in real space (as the density is increased) and the destruc- 
tion of the Fermi sea in momentum space (as the density 
is decreased). *I'hyb allows for an interpolation between 
momentum space {U* — 0) and real space (U* ^ 1). 
However, it never properly describes the Wigner crys- 
tal, since according to Bloch theorem, the 4>i{rj) are al- 
ways delocalized states (that can be concentrated around 
the UiS but that are delocalized anyway). The available 
values of momentum k are taken within the first Bril- 
louin zone, and hence, the liquid-hybrid transition can 
be viewed as an instability of the shape of the Fermi sur- 
face that goes from a circular to a hexagonal form. The 
symmetry is broken at this transition, but it is only in a 
second step that larger values of |fc| will come into play, 
allowing the (j)i {'fj) to get localized and the actual crystal- 
lization to take place. This transition in two steps, where 
first the direction of k and secondly its absolute value are 
affected, is reminiscent of the hexatic phase predicted in 
the classical melting. 

Stability of the hybrid phase. Fig^ shows the energy 
differences Enq — E^^y and i?hyb — £'cry as a function of 
Ts for a system of 72 electrons in a 48 x 84 grid. These 
energy differences are very small, less than 0.1% of the 
total energy of the sytem. rs ~ 40 where Enq — E^y « 
would be the critical value or Ts in the absence of the hy- 
brid phase 4^. However, we find that for 30 < Ts < 80, 
the hybrid phase has a smaller energy than both the liq- 
uid and the solid phase. Around r* « 30 we find a jump 
of U* from zero to U* — 0.3, see the inset of Fig. ^ 
U*/U « 0.015 up to rl* ^ 80, above which the crystal 
phase become more stable than the hybrid phase. Finite 
A'' corrections shown in the inset of Fig Q] at Ts w 42.2 do 
not perturb the previous picture. We note that although 
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FIG. 2: Energy of the liquid (circle, A = 6.0, B = 2.25), 
hybrid (triangle, A = 4.9, B = 2.58, U* = 0.3) and crystal 
(square, A — 4.9, B = 2.5, do — 2.95) phase as a function of 
imaginary time t for 72 electrons in a 48 x 84 grid at f/ = 20 
(rs ~ 42.2). Inset: variational energy of the hybrid phase as 
a function of (7* at (7 = 15 (r^ ~ 31.67 upper panel, A = 4.7, 
B = 2.3) and (7 = 30 (r, ^ 63.33 lower panel, A = 5.7, 
B = 2.9) 



the variational energy of the hybrid phase is lower than 
the one of the liquid it is still higher than the crystal 
variational energy. The FN-QMC treatment is necessary 
to show the stability of the hybrid phase, as shown in 
Fig. |2] To make contact with the calculations of |j, |y| , 
we have repeated these calculations for a more diluted 
system ly ~ 1/780 where the role of the underlying grid 
is negligible. The results are plotted in Fig. O in the 
same way as Fig. 2 of |^ (with ci = —2.2122). In the in- 
set of Fig^|3|we have reported the Wigner crystal data of 
[J| and jg for comparison and find a good quantitative 
agreement with the latter. We note that finite N cor- 
rections would raise the energies by an amount ^ 0.01 
(see the inset of Fig. P) and are difficult to evaluate, es- 
pecially in the absence of an analytical ansatz for the 
hybrid phase. 

Nature of the hybrid phase. It is important to real- 
ize that, as described above, the nature of the interme- 
diate phase is by construction something hybrid, being 
made of (delocalized) Bloch waves, yet having already 
the symmetry of the Wigner crystal. More insight can 
be gained by computing the electronic density (rif) (not 
shown) which is the superposition of peaks at the clas- 
sical positions (the Mi's) of the electrons in the crystal 
over a small background. The background is found to 
contain approximately 35% of the electrons while the 
rest lies in the peaks of the crystal. In that sense, the 
crystal part of the hybrid phase contains fewer electrons 
than sites, hence allowing exchange to take place. This 
is to be contrasted with the conjecture of Ref. [13 where 
a crystal with fewer sites than electrons was predicted. 
Although the total energy of the hybrid phase is be- 
low those of the liquid and crystal, both its kinetic and 
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FIG. 3; [2E{rs) — CiVs] jrj^ as a function of r^ for 72 elec- 
trons in a 180 x 312 grid. The curves for the liquid (circles), 
crystal (squares) and hybrid phase (triangle) can be compared 
directly with Fig. 2 of Q. Inset: same thing for the crystal 
phase only. The curves show our data for 72 electrons in a 
180 X 312 grid (squares) the data of Ref. 3 for 56 electrons 
(diamonds) and the data of Ref. Q for 56 electrons (trian- 
gles). 



electrostatic energies lie in between those of the liquid 
and solid. Fig. 0] shows the density-density correlation 
function (roughly measuring the probability of finding 
an electron at point r knowing that an electron is at 

ctctcjjcp^ for the three phases. 
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point 0), 5(7^) - ^(^_i) 
5hyb(^ for the hybrid phase is intermediate between a 
liquid and a crystal. The value of 5hyb(r) at its peaks is 
only twice as big as in the valley to be compared to a 
factor 15 at Ts — 100. In fact a very good fit is obtained 
with 5hyb(r) ~ 0.35 giiq(r) -I- 0.65 5cry(r)- 

To conclude, we find that a new quantum phase is to 
be expected instead of a direct melting of the Wigner 
crystal. This intermediate phase, whose physical proper- 
ties remain to be investigated in more depths, has hybrid 
behaviour between those of a solid and a liquid. 
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